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Abstract 

Background: Clear examples of ecological speciation exist, often involving divergence in trophic morphology. 
However, substantial variation also exists in how far the ecological speciation process proceeds, potentially linked to 
the number of ecological axes, traits, or genes subject to divergent selection. In addition, recent studies highlight 
how differentiation might occur between the sexes, rather than between populations. We examine variation in 
trophic morphology in two host-plant ecotypes of walking-stick insects {Timemo cristinoe), known to have diverged 
in morphological traits related to crypsis and predator avoidance, and to have reached an intermediate point in the 
ecological speciation process. Here we test how host plant use, sex, and rearing environment affect variation in 
trophic morphology in this species using traditional multivariate, novel kernel density based and Bayesian 
morphometric analyses. 

Results: Contrary to expectations, we find limited host-associated divergence in mandible shape. Instead, the main 
predictor of shape variation is sex, with secondary roles of population of origin and rearing environment. 

Conclusion: Our results show that trophic morphology does not strongly contribute to host-adapted ecotype 
divergence in T. cristinoe and that traits can respond to complex selection regimes by diverging along different 
intraspecific lines, thereby impeding progress toward speciation. 

Keywords: Sexual dimorphism, Timema cristinoe, Trophic morphology. Mandibles, Geometric morphometries, 
Bayesian clustering. Morphological uniqueness. Occupied morphospace. Disruptive selection. Selection dissipation 



Background 

Recent studies clearly demonstrate ecological causes for 
speciation [1-3], but extensive variance has also been 
noted in how far this process proceeds [4-6]. Speciation 
is often considered a continuum of divergence ranging 
from continuous variation among individuals within one 
species to complete discontinuity among groups in gen- 
etic, physiological and phenotypic adaptations to specific 
environments. The degree of divergence among taxa can 
be quantified along this continuum, typically by estimat- 
ing for example; levels of reproductive isolation, genetic 
dissimilarities or the amount of lineage sorting among 
groups [5-7]. This continuum analogy can lead to the 
misconception that all diverging taxa, given enough 

* Correspondence: denisroyl(S)gmail.com 

^Great Lakes Institute for Environmental Research, University of Windsor, 
Windsor, ON N9B 3P4, Canada 

^EAWAG, Swiss Federal Institute of Aquatic Science & Technology, 
Department of Fish Ecology & Evolution, Centre for Ecology, Evolution & 
Biogeochemistry, Seestrasse 79, CH-6047, Kastanienbaum, Switzerland 
Full list of author information is available at the end of the article 

(3 BioMed Central 



time, will eventually complete ecological speciation. 
However, divergence in many systems can reach points 
whereby some level of divergence has been achieved but 
is balanced by other factors that prevent it from moving 
further [5-8]. 

Observed variation in the progress towards ecological 
speciation may thus be due to differential divergence 
time, degree of gene flow, strength of divergent selec- 
tion, and the inherent genetic architectures underlying 
adaptive phenotypes [5-8]. However, another factor may 
be the number of ecological dimensions, traits or genes 
upon which divergent selection acts [3-6,8,9]. All else 
being equal, speciation may proceed further when diver- 
gent selection acts on a greater number of ecological 
axes, traits, or genes (multifarious selection, c.f. [4,5,9]). 
Selective pressures, however, can also act differently on 
different traits segregating them along different intraspe- 
cific factor (e.g., between sexes or between populations) 
or, alternatively, can lead to the expression of variable 
morphologies from a common genotype through the 
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development of phenotypic plasticity [10,11]. Thus, di- 
vergence can occurs to varying degrees among various 
traits, with some contributing more toward population 
divergence and speciation than others [4,5,8,11,12]. 

Of the possible traits involved in ecological speciation, 
feeding morphology has been shown to play a role in a 
wide array of taxa [2,13-18]. Theoretical work suggests 
that under conditions where the most common pheno- 
type depletes the most abundant resource, density- 
dependent disruptive selection on resource use traits can 
give rise to phenotypic variants adapted to the extremes 
of the resource spectrum [11,12,19,20]. This scenario 
favours extreme phenotypes that dissipate density- 
dependent selection regimes leading to resource 
partitioning among groups and, if these traits are also 
related to reproductive isolation, to resource based 
ecological speciation [11,12,19,20]. However, popula- 
tion divergence is not the only possible means of dissi- 
pating density-dependent selection. Density-dependent 
disruptive selection can also generate trophic traits 
that are sexually dimorphic, developmental stage spe- 
cific, or even highly plastic [10-12]. In such cases the 
development of reproductive isolation based on 
trophic trait divergence is less likely and may even be 
hindered from developing based on other traits as well 
[11,12,20,21]. Theory also posits that occurrence of 
any one of the above alternatives will alleviate selection 
and thus limit the impetus for further divergence in 
the same trait(s) along any other factors (e.g., sex, 
ontogenetic stage, population) [11,12,19-21]. Beyond 
theory, however, empirical studies are needed to clarify 
along which factors disruptive selection drives trophic 
morphological divergence. In particular, studies conducted 
in systems known to be undergoing ecological speciation 
can test whether trophic divergence contributes to the 
speciation process [4-6,22]. 

To address these questions concerning trophic diver- 
gence and ecological speciation we quantify mandible 
shape variation (known trophic traits in many insects; 
[16,23]) in multiple populations of T. cristinae that ex- 
hibit two ecotypes' found living on two different host- 
plant species. These ecotypes exhibit adaptive genetically 
based divergence in colour, colour-pattern, body size 
and shape [24-26] and also partial reproductive isolation 
[27-29]. Thus, we might expect host-associated diver- 
gence in trophic morphology to follow suit, although 
this has not been examined in past work. Here, we as- 
sess the relative contribution of four explanatory factors 
on mandible shape variation: (1) host-plant use in na- 
ture, (2) host rearing environment in the laboratory, (3) 
amount of gene flow into populations from populations 
on the alternate host, and (4) sex. Our design allows 
identification of both genetic and environmental (pheno- 
typic plasticity) effects on morphological variation. We 



assess shape variation in multi-dimensional morpho- 
logical space (morphospace) using traditional and more 
novel shape-based analyses. Contrary to expectations, 
we find little host-associated divergence in mandible 
shape. Instead, the main predictor of mandible shape 
variation is sex, with lesser roles of genetic background, 
rearing environment and host use. Although our results 
suggest that divergence in traits along one intraspecific 
factor may limit trait divergence along other factors, 
they also challenge the notion that divergence along one 
factor completely inhibits that expressed in others 
[12,19,20]. Our results thus provide new insights into 
ecological speciation in this system, and on the role of 
sexual dimorphism in the adaptive divergence of trophic 
morphology in general. 

Study system: Ecological speciation in Timema cristinae 

Timema cristinae is a species of wingless insect from the 
chaparral of southwestern North America whose individ- 
uals feed and mate exclusively on the host plants on 
which they live [30,31]. This species exhibits two eco- 
types defined by the host plant species on which they 
are found (i.e., 'Ceanothus and 'Adenostoma' ecotypes; 
see [29]). Previous work has shown heritable morpho- 
logical differences between ecotypes that have evolved in 
response to divergent selection for crypsis exerted by 
visual predators such as lizards and birds. In contrast to 
morphological adaptations for crypsis, divergent physio- 
logical adaptation to the different hosts is lacking, with 
both ecotypes exhibiting higher fecundity on the same 
host {Ceanothus) when predators are absent [29,32]. The 
ecotypes are in the process of ecological speciation, with 
pairs of populations on different host-plant species 
exhibiting stronger reproductive isolation than pairs on 
the same host [27,28,33]. Moreover the degree of both 
adaptive divergence and reproductive isolation is in- 
versely related to levels of gene flow between ecotypes 
[24,25,34]. Given this evidence for host-mediated eco- 
logical speciation and for constraining effects of gene 
flow on divergence, we might expect both processes to 
also influence trophic morphology. 

The two host plants used by T, cristinae offer very dif- 
ferent properties in general structure, physiology and 
leaf characteristics [35]. Ceanothus spinosus (hereafter 
C) is a tree-like Rhamnaceae with thick broad leaves, 
high water conductance and relatively low mechanical 
strength. Adenostoma fasciculatum (hereafter A) mean- 
while, is a bush-like Rosacea, characterized by thin tough 
needle-like leaves with low water conductance but high 
mechanical strength [35]. Given these differences in food 
resources, T, cristinae populations adapted to the differ- 
ent hosts might possess specific feeding morphologies 
adapted to these differences. This is especially likely con- 
sidering past work demonstrating mandible shapes in 
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other insects to be related to the toughness and hardness 
of their diets [23,36,37]. If this were the case in T. 
cristinaey it would demonstrate at least two different 
ecological dimensions (trophic ecology and predator 
avoidance) working in the same direction during eco- 
logical speciation. Alternatively, if trophic divergence is 
manifest primarily along non-host associated lines, then 
it will have diverged along a different trajectory relative 
to traits involved in crypsis and may stymie how far eco- 
logical speciation proceeds. 

Methods 

Sample collection & preparation 

A total of 200 sexually mature individuals were ran- 
domly chosen from previously collected and preserved 
samples of six populations (A^=25-50/population, Table 1 
for details) that are the focus of past and ongoing re- 
search [25,38], and which represent differential host use 
and levels of geographic isolation. Sexes can be unam- 
biguously differentiated in mature individuals by stark 
differences in external genitalia [30,39]. A subset of indi- 
viduals from each population (range = 3-25) were reared 
from first-instar to sexual maturity on either their native 
and/or alternate hosts under laboratory conditions, as 
described in previous work [25,38]. Left and right man- 
dibles were excised from specimens under a dissecting 
microscope, cleaned by bleach immersion and mounted 
onto scanning electron microscope (SEM) stubs in a 
specific orientation (Figure 1). Each individuals man- 
dible orientation was calibrated against an initial set. 



Table 1 Summary of Timema cristinae used in this study 



Population 


N total 


Host 


FLD 


LAB-^ 


LAB-C 


^Influence 








M 


F 


M F 


M 


F 




Lingual 


















HVA 


24 


A 


17 


2 




5 




15.97 


LA 


45 


A 


1 


3 


14 8 


11 


8 


0.21 


LRN 


18 


A 


11 


4 






3 


0.39 


OUTA 


24 


A 


16 


2 




5 


1 


119.89 


PRC 


26 


C 


15 


1 




10 




2.62 


VPC 


52 


C 


10 


2 


14 11 


12 


3 


19.46 


Occlusial 


















HVA 


23 


A 


17 


2 




4 




15.97 


LA 


46 


A 


2 


3 


13 8 


12 


8 


0.21 


LRN 


17 


A 


9 


5 






3 


0.39 


OUTA 


23 


A 


15 


2 




5 


1 


119.89 


PRC 


26 


C 


15 


1 




10 




2.62 


VPC 


49 


C 


10 


1 


14 11 


10 


3 


19.46 



Host: A-Adenostoma fasciculatum, C-Ceanothus spinosus; Rearing conditions: 
FLD = Field, LAB = Laboratory; ^- Calculations of influence are described in the 
text and equal a proxy for gene flow, with larger values reflecting more gene 
flow into a population from the alternate host. 



which we used as our experimental standard. Micro- 
adjustments to the placement of mandibles on the stub 
or to the ESEM stand were used to achieve standard 
orientation. Broken mandibles or ones whose final orien- 
tation was too different from the standard were removed 
from the analyses. To acquire better overall coverage of 
total shape, left and right mandibles were set in different 
perspectives accentuating different features. Left mandi- 
bles were placed in a lingual view emphasizing more 
functional features: distal incisor region (DIR), the prox- 
imal molar region (PMR) and structural support region 
(SSR; Figure lA). Right mandibles were placed in an 
occlusial view focussing on the SSR but also showing the 
incisor dents (ID; [36]; Figure IB). SEM images were 
generated using a Philips® FEI ESEM scanning electron 
microscope set at 60-70KV and a scale factor of 200 [im. 
Table 1 summarizes the final count of individuals used 
in each population under various rearing conditions. 

Geometric morphometries 

Electron micrographs were subjected to geometric mor- 
phometric analyses using 11 homologous landmarks 
chosen to outline the features listed above. Left (here- 
after lingual) and right (hereafter occlusial) mandibles 
were analysed separately. Standardized landmark config- 
urations were imported into MorphoJ [40], where a 
generalized least squares Procrustes superimposition 
(GPS) was performed. The GPS generated a consensus 
configuration for both mandibles and a thin-plate spline 
was used to quantify deviation of individual specimen 
landmark configurations from the consensus along prin- 
cipal axes of shape change [40,41]. Deviations from the 
consensus along principal axes generate partial warps 
quantifying individual based shape changes [41,42]. For 
both mandibles, the thin-plate spline produced 16 partial 
warps and two uniform components per individual, 
representing local and overall shape variables, respectively. 
All 18 shape variables were used in subsequent analyses. 

Shape variable size correction 

Centroid sizes of the lingual mandibles were transformed 
to their natural logarithm to fit a normal distribution while 
no such transformation was required for the occlusial. 
Centroid sizes were tested for differences among four ex- 
planatory factors: sex, host-adapted ecotype, population of 
origin, and rearing condition. In both cases, females 
exhibited significantly larger mandibles than males (lingual; 
females = 992.67±64.84(im vs. males = 876.75 ± 50.83(im; 
^68.66 = -11.13, P < 0.001: occlusial; females = 953.35 ± 
67.83(im vs. males = 850.87 ± 44.76(im; ^9.91 = -9.56, P< 
0.001). No significant differences were found in centroid 
sizes with respect to other tested factors. Because size is 
known to influence shape through allometric trajectories 
[41,43,44], shape variables were regressed against centroid 
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Figure 1 Electron micrographs of Timema cristinae mandibles showing features outlined by landmarks used. (A) Lingual view sliowing 
distal incisor region (DIR), the proximal molar region (PMR) and the structural support region (SSR). (B) Occlusial view showing the incisor dents 
and the SSR. Consensus configuration of all specimen corrected landmark configurations of the lingual and occlusial views are superimposed on 
mandible figures. 



size using MANCOVA approach and residuals were used 
in subsequent analyses [40,44,45]. Multivariate size correc- 
tion used a single regression for both sexes because the 
interaction term between sex and size (used as covariate) 
did not account for significant shape variation {P > 0.05). 

Bayesian shape analyses and morphological clustering 

Using all data for each mandible orientation (Ungual 
N= 189; occlusial N= 184), we assessed the most likely 
number of statistically different groups within samples 
based solely on their shapes using a Bayesian posterior 
probability assignment tests run in the program 
Autoclass-C V 3.3.4. This analysis allowed us to statistically 
find the most obvious separation in the data without prior 
grouping by assessed factors. Detailed descriptions of 
Autoclass-C are available from the NASA home page 
(http://ti.arc.nasa.gov/tech/rse/synthesis-projects-appU- 
cations/autoclass/autoclass-c/) and in Cheeseman & 
Stutz (1996). Briefly, Autoclass-C uses a Bayesian ex- 
tension of finite-mixture modelling to perform un- 
supervised searches recovering the most probable 
number of statistically different groups given the data. 
Searches make no prior assumptions of actual group 
number but assume each resulting group fits a given 
distribution set by the user [46]. In each iteration, 
Autoclass-C generates a number of hypothetical distri- 
butions with given parameters (e.g., means and va- 
riances) to which actual data are permutated and fit 
over a given number of cycles. Convergence is achieved 
in each iteration when the actual data fit the hypothe- 
sized distributions within a given error estimate over a 



predefined number of cycles. The probability of the con- 
verged data is then evaluated using a Bayesian framework. 
We used 10 000 000 iterations allowing each to reach con- 
vergence over 100 000 cycles where convergence was 
deemed acceptable when actual and hypothesized parame- 
ters were within 0.0025 over at least 10 consecutive cycles. 
Searches recorded the most probable number of groups in 
the data every 10 000 iterations and saved the best 200 
overall. We assumed variables used in the modelling fit 
normal or lognormal distributions and that variable- 
specific error terms were fixed. Variable-specific error 
terms were calculated using all individuals included in the 
analyses. Autoclass-C was initially run using all shape vari- 
ables but was also run using the same settings on relative 
warps (RWs) (PCs of unweighted shape variables; [41]). Re- 
sults from both analyses were identical and so those using 
RWs, which fit model assumptions better, are presented. 
Autoclass-C results also generated individual-based poster- 
ior probabilities of belonging to recovered groups which 
were then used to generate probability of assignment plots 
[47]. Mean RW scores of individuals belonging to domin- 
ant clusters were also used to generate deformation grids 
outlining group specific shapes. 

Uniqueness of occupied morphospace 

Multivariate parametric analyses are most reliable when 
sample sizes among and within grouping factors are well 
balanced [41,48]. Such balanced designs, however, can 
be difficult to achieve for small complex morphological 
features (~ 620 [im total length), which are delicate, 
costly to prepare and not easily replaced. Because our 
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overall data set was not conducive to parametric assess- 
ments of variance partitioning, we chose to compare 
shape differences within explanatory factors by assessing 
morphological uniqueness. Morphological uniqueness 
(hereafter MU) quantifies the amount of unique mor- 
phological shape space occupied by two predefined 
groups. MU is based on the non-parametric niche 
overlap index developed to estimate the overlap between 
groups based on quantitative functional traits [49]. 
Briefly, along each RW, each individuals score is 
converted to a kernel distribution which contributes to 
an overall kernel density function formulated for the 
group to which it belongs [49]. Group specific functions 
for each RW are then compared by stepwise integration 
of the intersecting area between the two functions over 
the predefined range given by the maximum range of 
the largest group. This integral determines the overlap 
between the two groups along this particular RW 
[49,50]. Because the functions are bounded over the 
same range, the uniqueness along a RW can be consid- 
ered unity minus the overlap. The uniqueness calculated 
over each RW is then weighted by the amount of vari- 
ance accounted for by each RW (determined from 



eigenvalues) and summed. Generated MU indices range 
between 0 and 1 quantifying non-overlapping morphospace 
occupied by the two compared groups. Significance of MUs 
between groups is assessed by resampling the data for the 
same number of individuals but arbitrarily assigned to the 
different groups (with replacement) one thousand times. 
Because the MU can be applied to the same data but test- 
ing different explanatory factors, it can quantify how each 
explanatory factor partitions the same morphospace de- 
fined by a given variable (here RWs) or set of variables 
outlining more general morphospace. Here, overall MUs 
were calculated and compared among all groups within an 
explanatory factor (i.e., sex, host, population, and rearing 
condition) in a pairwise fashion and permutations were 
performed using R scripts (available from the authors; [51]). 
MUs for each explanatory factor were also calculated for 
each RW to assess which factor best partitioned the vari- 
ance along each. In analyses including factors with more 
than two groups (i.e., population and rearing condition), 
the mean MUs from all pairwise comparisons along each 
RW were used in comparisons with other factors (see 
Table 2). Deformation of the consensus configuration along 
the first three RWs was also determined by regressing 



Table 2 Decomposed raw Morphological Uniqueness of T. cristinae mandibles quantified among all tested factors 



RW 


%a 


Sex 


Host 


Population 


Rearing 
condition 


Dominant 
factor 


%a 


Sex 


Host 


Population 


Rearing 
condition 


Dominant 
factor 




Lingual 












Occlusiol 












1 


31.38 


0.087 


0.164 


0.299 


0.201 


P 


27.81 


0.514 


0.177 


0.310 


0.231 


S 


2 


14.62 


0.511 


0.149 


0.255 


0.301 


S 


14.25 


0.112 


0.105 


0.290 


0.216 


P 


3 


11.02 


0.156 


0.300 


0.259 


0.218 


H 


13.90 


0.180 


0.192 


0.196 


0.164 


P 


4 


8.92 


0.565 


0.134 


0.231 


0.200 


S 


8.19 


0.410 


0.069 


0.274 


0.140 


S 


5 


7.60 


0.195 


0.098 


0.218 


0.140 


P 


6.86 


0.124 


0.088 


0.151 


0.212 


R 


6 


6.04 


0.757 


0.180 


0.274 


0.232 


S 


6.05 


0.167 


0.121 


0.240 


0.164 


P 


7 


4.91 


0.463 


0.127 


0.239 


0.143 


s 


4.72 


0.522 


0.066 


0.283 


0.250 


S 


8 


4.19 


0.076 


0.114 


0.151 


0.179 


R 


4.28 


0.491 


0.144 


0.194 


0.170 


S 


9 


3.13 


0.675 


0.150 


0.226 


0.112 


S 


3.21 


0.150 


0.120 


0.217 


0.101 


P 


10 


2.01 


0.503 


0.148 


0.250 


0.201 


S 


2.86 


0.519 


0.109 


0.173 


0.136 


S 


11 


1.42 


0.383 


0.082 


0.196 


0.129 


S 


2.33 


0.212 


0.074 


0.165 


0.217 


R 


12 


1.23 


0.117 


0.113 


0.196 


0.195 


P 


1.47 


0.677 


0.118 


0.195 


0.191 


S 


13 


1.00 


0.130 


0.106 


0.215 


0.120 


P 


1.08 


0.672 


0.148 


0.238 


0.227 


S 


14 


0.71 


0.226 


0.074 


0.216 


0.183 


S 


0.86 


0.362 


0.065 


0.193 


0.149 


S 


15 


0.64 


0.374 


0.175 


0.254 


0.201 


S 


0.67 


0.380 


0.115 


0.218 


0.129 


S 


16 


047 


0.203 


0.210 


0.206 


0.167 


H 


0.61 


0.171 


0.109 


0.216 


0.167 


P 


17 


0.39 


0.136 


0.071 


0.182 


0.166 


P 


0.50 


0.142 


0.078 


0.189 


0.144 


P 


18 


0.30 


0.129 


0.125 


0.199 


0.170 


R 


0.33 


0.315 


0.153 


0.256 


0.110 


s 






9 


2 


5 


2 






10 


0 


6 


2 





Row entries describe divergence along each Relative Warp (PC of shape variables) determined from pairwise comparisons within each tested factor. Population 
and rearing condition entries show means calculated over 15 and 6 comparisons, respectively. Calculated using all individuals (A/= 189 lingual and A/= 184 
occlusial). H = host, P = population, R = Rearing condition, RW = Relative warp, RWdom - number of dominant relative warps (does not consider the variance 
explained by specific RWs), S = sex,% o = percent variance explained along RW. 
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individual shape variables onto RWl, RW2 and RW3 
scores, respectively [40,52,53]. 

Genetic and environmental basis of mandible shape 

To quantify the relative importance of genetic back- 
ground (e.g., population of origin) and rearing environ- 
ment on mandible shape (s), we used individuals from 
populations LA and VPC reared from first instar to sex- 
ual maturity on both hosts in a reciprocal transplant ex- 
periment performed in the laboratory (Table 1; [25]). 
Landmark configurations for these samples were isolated 
and used to generate a new consensus configuration and 
a new set of size corrected shape variable as described 
above. We then used MANOVA to determine the influ- 
ence of sex, population of origin, host rearing condition 
and their interactions on mandible shape variation. The 
population effect tested influences of genetic background 
while host rearing condition tested environmental 
effects. Interaction terms tested sex specific responses in 
the different populations and hosts, and population 
specific responses to hosts [25,54]. MANOVA was 
performed using default settings and a type III sum of 
squares model taking into account uneven sample sizes 
in STATISTICA 64 v.lO (Stat Soft Inc., Tulsa, OK USA). 

Sex removal 

Results of the analyses revealed an overwhelmingly 
strong effect of sex on mandible shape (see results). 
However, other factors also explained some portion of 
mandible shape variation. To justify the consideration of 
a single sex (male) in subsequent analyses and to relate 
findings to the entire population, we compared the 
shape changes occurring in all males and females in both 
common and global morphospace using ordered axes 
analysis [55] and the assessment of angular differences 
along RWs between sexes [41,55]. Results from both 
analyses indicated that males and females, although 
differing in overall shape, followed the same shape 
change trajectories through both common and global 
morphospace (additional file 1). Consequently, data were 
sub-divided by sex and Bayesian assignments, MU and 
MANOVA tests described above were performed using 
only field-collected (wild) males from populations with 
nine or more individuals. In this case, MANOVAs 
nested populations within host as specific populations 
use only a single host in the wild. 

Gene flow 

Finally, we tested the influence of gene flow between 
hosts-adapted populations on mandible shape variation. 
Isolated populations with low levels of gene flow might 
be expected to exhibit less shape variation than those 
characterized by higher gene flow. Using only wild males 
from populations with nine or more individuals (Table 1), 



we generated another consensus configuration and set of 
shape variables. Shape variables were converted to RWs 
as above and population-specific standard error of the 
mean of RWl -RWS scores (an index of shape disper- 
sion) were regressed against a population-specific index 
of influence, a proxy for gene flow. The index of influ- 
ence for each population was calculated as the sum of 
the proportion of the area of the focal population and 
its adjacent population occupied by the alternate host 
(determined from aerial photographs; [56]) divided by 
the sum of the distance between them. Previous 
work has demonstrated that the proportion of aerial 
coverage occupied by the alternate host is highly corre- 
lated with levels of gene flow inferred from molecular 
data [25,28,56,57]. 

Results 

Bayesian assignment by shape 

From the overall dataset of both mandibles, the most 
likely number of independent groups recovered from 
unsupervised Bayesian searches separated individuals 
into three clusters with two dominating the assignments 
and one encompassing only three individuals (Figure 2). 
In both cases, the two dominant clusters were strongly 
associated with sex indicating that sexes were the most 
obvious separation in the data. The third minor cluster 
was likely artifactual resulting from Autoclass-Cs sensi- 
tivity in generating groups with such small fixed error 
terms, as its membership was not consistent between 
mandibles. In lingual mandibles, 97% of males were 
assigned to one cluster with greater than 95% posterior 
probability, while 96% of females were assigned to a dif- 
ferent cluster also with greater than 95% posterior prob- 
ability (Figure 2A). The best classification was 2.64 times 
more probable than the next best classification scheme 
separating individuals into the same clusters and 104.27 
times more probable than the classification estimating 
the next most probable number of clusters (two; male 
and female only). Assignments for occlusial mandibles 
were similar to those of the lingual, but varied in cluster 
statistics (Figure 2B). Ninety-one percent of females 
were assigned to the same cluster with greater than 95% 
posterior probability. Ninety-eight percent of males were 
assigned to a single cluster with greater than 95% poster- 
ior probability. The final best classification was 163.04 
times better than the next most probable classification 
and 160 000 times better than the classification estimat- 
ing the next most probable number of clusters (two; 
male and female only). Lingual deformation grids 
showed that for the same overall sized mandible, the fe- 
males was generally narrower with a smaller proximal 
molar region (PMR), structural support region (SSR) and 
more pointed distal incisor region (DIR) (Figure 2A). 
Average male occlusial mandibles appeared to have more 
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females 




males 
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Figure 2 Bayesian posterior probability plots estimating number of clusters among sampled individuals based on lingual (A) and 
occlusial (B) mandibles shapes. Each individual is represented by a bar wliose colouration is based on its probability of belonging to each 
cluster. Deformation grids associated with the two main clusters derived from the average shape of individuals belonging to each. 



streamlined SSR and longer more extended incisor dents 
(ID) (Figure 2B). 

Factor-specific MU on overall data 

In both mandibles, sex best partitioned the overall 
morphospace assessed along all 18 RWs (overall MU), 
followed by mean pairwise differences between popula- 
tions and rearing conditions, while host-adapted eco- 
types partitioned the morphospace the least (Figures 3, 
Additional file 1: S2 and S3). Although sex best 
partitioned the morphospace along a majority of individ- 
ual RWs (RW-specific MU), other factors also 
partitioned morphospace along the same RWs, but to a 
lesser degree (Table 2). This result supports the notion 
that all factors interact to some extent but that one 
typically dominate the divergence (Table 2). After sex, 
mean pairwise population comparisons best partitioned 
morphospace along a larger number of RWs followed by 
mean pairwise rearing condition comparisons and then 
host-adapted ecotypes (Table 2). In the morphospace 
outlined by the first three RWs, sexes were clearly most 
separated along RW2 for the lingual mandible but along 
RWl for the occlusial. In both cases shape differences 
along these RWs, as depicted from deformation grids 
along RWs, followed similar trends as that determined 
in the Bayesian clustering (compare Figures 2 and 3). 
When populations and rearing conditions occupied 
significantly different morphospace, separation was typ- 
ically most obvious along RWl and RW2. However, 
many populations also showed substantial overlap con- 
sistent with non-significant MUs (e.g.. Figure 3C,F and 
Additional file 1: Figures S2, S3). In both mandibles, 
populations sharing the same host tended to have larger 



MUs than those on alternate ones (MU same = 27, 
different = 24.75). Similarly, the largest observed differ- 
ences in MU among rearing conditions occurred be- 
tween those both using A as host (i.e., lab versus field 
environments), while that quantified on C were not sig- 
nificant. Thus, in both population and rearing condition 
comparisons, anticipated higher MUs on alternate hosts 
was not observed, indicating a low level of host influence 
in these factors. Despite their significant MUs, no obvi- 
ous separation was observed along any of the three RWs 
when considering host-adapted ecotypes (Figure 3B,E), 
indicating limited host-specific shape adaptation when 
all data are considered together. 

Genetic and environmental basis of mandible shape 

MANOVA results from the reciprocal transplant experi- 
ment revealed a pattern similar to the overall dataset. 
Namely, sex accounted for most of the shape variation 
for both mandible orientations (Table 3). However, sig- 
nificant population (genetic) and host rearing condition 
(environmental) components were also detected. For the 
lingual mandible, genetic effects accounted for a greater 
portion of the variance than did environmental ones, 
while for the occlusial, both factors accounted for similar 
variance proportions. Both mandibles also demonstrated 
a significant interaction in shape variation between 
population and host implying that one population was 
more responsive to environmental differences than the 
other (Table 3). No sex specific influences of population 
or host rearing condition were observed in either mandi- 
bles indicating that both sexes responded similarly to 
genetic and environmental differences. The power of 
these interactions to reject the null, although low, was 
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Figure 3 Morphological uniqueness (MU) assessed for T. cristinae lingual and occlusial mandibles using all 18 relative warps (RWs) 
supported by 1000 resampling permutations. Panels show 80% confidence bubbles outlining morphospace occupied by different A) sexes B) host 
plant ecotypes and C) rearing conditions for the lingual mandible along the first three RWs accounting for 57% of overall shape variation. Panels D-F show 
the same as A-B but for the occlusial mandible whose first 3 RWs account for 56% of the shape variation. Deformation grids next to RWs in panel A and F 
show general shape change trends along each. FLD = field LAB = Laboratory A = Adenostoma and C= Ceonothus; §-Mean MU calculated from all pairwise 
comparisons; *** = P< 0.001 ,** = P< 0.01 , and * = P < 0.05. 



nevertheless above 35% in all cases [48], and results are 
consistent with those of the ordered axes analyses 
showing similar shape change trajectories in both sexes 
regardless of other factors in the overall data (see 
Additional file 1: Figure SI). 

Host- and population-specific divergence within sex 

Bayesian classification schemes of field-collected (wild) 
males based on shape variables for both mandibles re- 
covered only a single most probable cluster. For the lin- 
gual and the occlusial mandibles, best classifications 
were 8 170 000 and 10.15 times more likely than the 
next best classifications predicting two cluster, respect- 
ively. Thus, Bayesian analyses could not further cluster 
the data when only a single sex was considered. In con- 
trast to Bayesian results, MANOVA and MU analyses 
showed variable effects of host and population on 
mandible shapes. For the lingual mandible, MANOVA 
showed that host accounted for a larger portion of the 
overall variance than did populations nested within host 
(Table 4). For the occlusial mandible, though, host was 
not a significant explanatory factor, but population 



within host was (Table 4). MU estimated along all 18 
RW did not agree with MANOVA results but rather 
showed that both host and population significantly 
partitioned outlined morphospace for both mandibles. 
In both cases, however, the MU calculated between 
host-adapted ecotypes was less than that expressed be- 
tween individuals belonging to different populations, es- 
pecially when only significant population comparisons 
were considered (Figure 4 and Additional file 1: Figures 
S4-S5). For the lingual mandible, although the MU was 
significant using all 18 RWs, no clear separation between 
host-adapted ecotypes was obvious in the morphospace 
outlined by the first three RWs, indicating a lack of 
host-specific shape as defined by the deformation grids 
along RWl-3 (Figure 4A). Population level MUs, mean- 
while, were stronger and more obvious along RW2 
where some populations had more pointed DIRs, more 
curved PMRs and broader SRRs than others (Figure 4B-C 
and Additional file 1: Figure S4). Half the significant popu- 
lation MUs recovered were between populations sharing 
the same host indicating that not all population differ- 
ences were host-based. For the occlusial mandible, clear 
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Table 3 MANOVA testing the genetic and environmental basis of mandible shape in Timema cristinae (using LA and 
VPC reared under laboratory conditions) 



Sources of variation 


Wilk's 
A 


F 


df 

Ho/ Error 


P 


rciilicii Vaiiaii^c 

explained 

(% n') 


Observed 
power 


Lingual 
















Sex 


0.19 


13.15 


1 9. 
1 o 


J / 


<0.001 


80.6 


1.00 


Population (Popn) 


0.25 


9.30 


1 8 
1 O 


J / 


<0.001 


74.6 


1.00 


Host rearing condition (HRC) 


0.54 


2.65 


1 8 


57 


0.003 


45.6 


0.99 


Sex X Popn 


0.65 


1.68 


18 


57 


0.072 


- 


0.87 


Sex X HRC 


0.82 


0.70 


18 


57 


0.794 


- 


0.42 


Popn X HRC 


0.49 


3.26 


18 


57 


<0.001 


50.7 


0.99 


Occlusial 
















Sex 


0.28 


7.69 


18 


55 


<0.001 


71.6 


1.00 


Popn 


0.36 


5.49 


18 


55 


<0.001 


64.3 


1.00 


HRC 


0.35 


5.68 


18 


55 


<0.001 


65.0 


1.00 


Sex X Popn 


0.78 


0.85 


18 


55 


0.682 




0.36 


Sex X HRC 


0.84 


0.60 


18 


55 


0.855 




0.52 


Popn X HRC 


0.44 


3.83 


18 


55 


<0.001 


55.6 


1.00 


q^- values indicate tlie amount of variance attributable to tested factors; 

F - approximation to F-ration derived from reported Wilk's A and degrees of freedom. 



differences were observed along RWl between host- 
adapted ecotypes where A-ecotypes scored more positively 
and exhibited more slender SSRs and longer IDs relative 
to C-ecotypes (Figure 4D). Populations exhibiting signifi- 
cant MUs were also most obviously different along RWl, 
where populations using A also scored more positively 
and shared the same features as A -ecotypes (Figure 4E-F 
and Additional file 1: Figure S5). All significant population 
MUs in the occlusial mandible were between those using 
alternate hosts suggesting a strong host-population inter- 
action (Figure 4E-F and Additional file 1: Figure S5). 

In summary, single sex analyses revealed variable levels 
of population and host mandible shape differences de- 
pending on the test performed (Bayesian, MANOVA or 
MU) likely reflecting more subtle shape differences than 
those observed between sexes. Pairwise population 
MUs were generally greater than those between host- 



adapted ecotypes, but host-population interactions, es- 
pecially in the occlusial mandible, were recovered. 
Thus, host effects on mandible shape were slight and 
subtle at best. 

Gene flow 

Mandible shape dispersion estimated from wild males 
was not significantly related to the potential for gene 
flow among neighbouring populations for either man- 
dible orientation {P > 0.05; Additional file 1: Figure S6). 

Discussion 

Here, we showed mandible shape differentiation between 
sexes and between host- adapted ecotypes, some of which 
was population specific. More slender occlusial mandi- 
bles with longer more extended cutting regions were 
associated with individuals using Adenostoma as host. 



Table 4 MANOVA testing mandibles shape variables among field collected male T. cristinae 



Sources of variation 



Wilk's 
A 



df 



Ho/ Error 



Partial variance 
explained {%, r^) 



Observed 
power 



Lingual 
















Host 


0.34 


4.96 


18 


47 


<0.001 


65.5 


0.99 


Population (Host) 


0.18 


2.02 


54 


140.86 


<0.001 


44.4 


0.99 


Occlusial 
















Host 


0.63 


1.38 


18 


43 


0.189 




0.76 


Population (Host) 


0.18 


1.90 


54 


128.94 


0.002 


44.1 


0.99 



rj - values indicate the amount of variance attributable to tested factors; 

F - approximation to F-ration derived from reported Wilk's A and degrees of freedom. 
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Figure 4 ML) of field reared male T. cristinae lingual and occlusial mandible shapes assessed using all 18 RWs supported by 1000 
resampling permutations. Panels show 80% confidence bubbles outlining morphospace occupied by different A) host plant ecotypes, B and C) key 
pairwise populations comparisons along the first three RWs accounting for 60.3% of shape variation. Panels D-F show the same as A-C, but for the 
occlusial mandible whose first 3 RWs account for 67% of the shape variation. Deformation grids next to RWs in panel C and D show general shape 
change trends along each. = P < 0.001 , = P < 0.01 , and * = P < 0.05. 



whose resources are likely more difficult to access [35]. 
However, we also showed that the strongest predictor of 
mandible shape variation in this species is sex rather 
than host-adapted ecotype, population of origin or rear- 
ing condition. This sexual dimorphism was prevalent both 
in the overall data and in laboratory-reared (reciprocal 
transplant) individuals from both hosts, indicating its 
heritable nature. The degree to which genetics and the en- 
vironment influenced mandible shape was variable, and 
depended on the mandible features emphasized in the two 
mandible orientations (lingual vs. occlusial). Once sexual 
dimorphism and rearing conditions were accounted for 
(wild male data), we showed that host played a significant 
yet variable role in determining mandible shapes, as did 
population of origin. Finally, and in contrast to most other 
morphological features measured in this species to date, no 
effect of gene flow was observed on mandible shape vari- 
ation. Below, we further develop these issues and discuss 
their evolutionary and ecological implications. 

Potential causes of sexual dimorphism in T. cristinae 
mandibles 

Sexual dimorphism is common among various animal 
taxa and is often attributable to divergent sexual and/or 
ecological roles between the sexes [58-61]. The causes of 
sexual dimorphism in trophic morphology in T. cristinae 



remain unknown but could be attributable to various 
mechanisms, including fecundity-based sexual size di- 
morphism, sexual selection, or differential feeding ecol- 
ogy between sexes. 

Sexual size dimorphism 

Sexual size dimorphism can be a consequence of adapta- 
tion to different male and female reproductive roles 
[58-60]. These roles are often the result of fecundity se- 
lection on egg development, laying and nutrient alloca- 
tion amongst others. For instance, female Timema, 
unlike males, ingest soil from below their host plant, 
which they use to coat egg cases. Consequently, the egg 
laying sex often grows faster and larger to compensate 
for the extra allocation of resources [59,60]. Because 
such traits are costly to express in the absence of fe- 
cundity selection, such fecundity-based adaptation can 
translate into size differences between sexes and into di- 
morphism based solely on allometry [59,60]. Previous 
work has demonstrated sexual dimorphism in other 
morphological traits of T, cristinae in relation to size, 
with female body size being larger than males [25,39]. 
Here, we found female mandibles to be significantly lar- 
ger than those of males, consistent with sexual size di- 
morphism. However, strong sexual dimorphism was stiU 
present in the residuals of shape variable regression with 
centroid size, a procedure that should account for 
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allometric differences between sexes [40,44]. Thus, the 
sexual dimorphism observed here goes beyond simple 
allometric differences. Moreover, in most reported in- 
stances of fecundity-based sexual size dimorphism, sexu- 
ally dimorphic traits consistently co-vary with size 
differences (i.e., the larger sex also exhibits larger traits) 
[58,62]. Instead, we found that size-corrected mandibles 
showed much broader (lingual) and more extended 
(occlusial) features in males than in females, opposite to 
the overall dimorphic body size-trait covariance pattern. 
Thus, other factors than size are involved in the sexually 
dimorphic mandible shapes in this species. 

Sexual selection 

Sexually dimorphic mandibles may be the result of sex- 
ual selection in efforts to sequester mating opportunities 
from conspecific rivals [59,63,64]. Within insect taxa in 
particular, male fighting behaviours can involve the use 
of mandibles [63,65]. Reports of male fighting in 
Phasmatodea (including T. cristinae) however, are rare 
but can arise during mate guarding where a male en- 
counters an already coupled pair [66,67]. In these en- 
counters, coupled males grip the female abdomen with 
their genitalia curving it away from approaching rivals to 
prevent copulation. Males may also engage in a 'boxing' 
behaviour using their front legs [67], but the use of man- 
dibles in such interactions has not been reported. In 
addition, more recent T, cristinae mating behaviour 
work does not show any evidence of mandible use in 
mate pairing, courtship display, initiation of copulation, 
or post-copulatory mate guarding [68]. Thus, although 
the usefulness of mandibles as sexual weapons or traits 
conferring reproductive success in T, cristinae has not 
been explicitly tested, sexually dimorphic mandibles do 
not appear to be related to sexual selection in any obvi- 
ous way. 

Ecologically based sexual dimorphism 

Another explanation for sexually dimorphic mandibles is 
the development of sex-specific feeding ecologies, as 
reported in other species [69-72]. Although T, cristinae 
males and females carry out all life history stages on the 
same plant species, it may be that they partition food re- 
sources by feeding on different parts of the same host- 
plant. Such tissue specific feeding behaviours have been 
demonstrated among other plant-browsing insects 
[16,73], but not typically between sexes within the same 
species. We showed that female lingual mandibles were 
more slender and pointed relative to those of males. 
Studies have demonstrated that sharper more pointed 
mandibles are generally better at initiating and propagat- 
ing fractures in tough leaf material [16,73]. Thus, female 
mandibles may be better adapted to feeding on the 
tougher parts of host plants whereas male mandibles 



may be better on the softer mesenchyme leaf tissue. 
Such sex-specific feeding strategies are also consistent 
with observed differences in occlusial mandibles wherein 
females had broader structural features with more 
curved cutting regions than males. Females may house 
larger mandibular muscles in enlarged structural regions 
capable of applying greater pressure to sharper cutting 
features enabling easier plant tissue fracture. Sex-specific 
feeding habits have not been tested in T, cristinae^ but 
offer clear avenues for further research. Directed studies 
quantifying sex-specific diets using either dietary tracer 
information (i.e., stable isotopes and/or fatty acids) or 
tissues enriched with other chemical tracers could help 
determine whether mandible sexual dimorphism is diet 
based [74,75]. 

Sex versus host and other factors 

The large sexual dimorphism in T, cristinae mandible 
shape suggests that it may limit mandible divergence 
along other intraspecific factors such as host or popula- 
tion. Theoretical modelling suggests that the mechanism 
by which adaptive ecologically-based sexual dimorphism 
emerges may be the same as that driving adaptive eco- 
logical speciation or phenotypic plasticity [10-12,19,20]. 
In this context, the development of ecologically- 
based dimorphic sexes, ecologically divergent species, or 
ecology-based adaptive phenotypic plasticity, is thought 
to supplant the evolution of other alternatives by dissi- 
pating disruptive selection among traits [11,12,19,20]. 
The development of any one of these types of divergence 
should thereby alleviate the impetus for other types of 
divergence in the same trait(s) [10-12]. If mandible 
shape differences between males and females result from 
adaptation to different diets, then the mandible shape 
adaptations between sexes might limit further host- 
based (or other factor based) mandible shape adaptation 
(but see [12]). The large sexual divergence in mandible 
shapes observed in both the overall data and in the re- 
ciprocal transplant experiment attests to the dominance 
of sex in determining mandible shape. However, the lack 
of sex-specific effects of host rearing condition or popu- 
lation in the reciprocal transplant experiment indicate 
that genetic and environmental factors tend to influence 
both sexes similarly (i.e., that T, cristinae do not exhibit 
sex-specific adaptations to the different host plants and 
that the sexes do not feed on different hosts). This is fur- 
ther supported by the fact that although males and fe- 
males accumulate variance in mandible morphology at 
different evolutionary rates, they do so along the same 
shape change trajectories in both common and global 
morphospace (see Additional file 1). Thus, mandible 
shape differences between sexes are not host plant- 
specific. Whether the sexes feed on different parts of the 
same host plant, however, is unclear (see above). 
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Moreover, and although sex best partitioned the vari- 
ation in mandible shape, it did not totally preclude man- 
dible shape divergence along other intraspecific factors 
along the same axes of shape variance (RWs in Table 2), 
a feature inconsistent with the theoretical treatments 
above (but see [21]). This may be because the theoretical 
models described assume competition for resources gener- 
ate the disruptive selection that drives trait divergence. It 
remains unclear, however, whether such a scenario applies 
to T. cristinae, as food resources seem readily available. 
Nevertheless, and regardless of the cause, mandible shape 
divergence between the sexes dominates that observed be- 
tween host-adapted ecotypes and/or among populations. 
These findings contrast previous work in this species 
reporting strong divergence between host-adapted eco- 
types in most other morphological features [4,25,27,29,38]. 
Thus, our data suggest that trophic morphology (as quanti- 
fied here) does not strongly respond to the multifarious se- 
lection regime driving host-based ecotype adaptation, nor 
contribute to reproductive isolation between host-adapted 
ecotypes. It may instead limit this process by occurring 
along a very different intraspecific factor. 

Host and population based shape In wild males 

Once the influence of sex was removed, mandible shape in 
wild males along a gene flow gradient was variably explained 
by different factors depending on the test performed. Dis- 
crepancies are likely the consequence of differential test sen- 
sitivities and the more subtle nature of within sex shape 
differences. MANOVA used all RWs weighted equally to es- 
tablish groups differences, and may therefore be oversensi- 
tive to small differences in less important RWs. Likewise, 
Bayesian clustering also considered all RWs evenly (al- 
though see [76]) and found the most likely clustering based 
on normal or log-normal distributions with fixed error 
terms. If differences between groups are more subtle and 
organised in complex hierarchies, they may be more difficult 
to recover [46,76]. In contrast, the kernel density based MU 
weighs the contribution of each RW to group-specific differ- 
ences in occupied morphospace, and may thus more accur- 
ately quantify morphological differences [49]. MU indices, 
however, cannot easily quantify interactions among factors 
and can therefore underestimate how they influence overall 
shape [48]. Nevertheless, and taking into account these 
discrepancies, lingual mandible shape variation was not con- 
sistent with host-based adaptive feeding. If it were, A-eco- 
types would be expected to develop shapes that more easily 
initiate and/or propagate fractures in tougher more resistant 
materials [16,35,73]. Pointed sharp blades would do this best 
[16,73]. Contrary to these expectations, and although MU 
results showed some morphospace partitioning between 
host-adapted ecotypes, this was not well reflected in shape 
differences outlined by the first three RWs. Moreover, pair- 
wise population level comparisons recovered as many 



significant MUs between populations using the same host as 
those using alternate ones. Thus lingual mandible shapes, 
chosen to reflect more functional features, were likely least 
influenced by host level adaptations. Variation in the 
occlusial mandible, however, showed significant levels of 
host-adapted ecotype variation. Separation in occupied 
morphospace by the different host-adapted ecotypes was 
clear along RWl, where A-ecotypes had shapes more appro- 
priate for a harder, tougher host plant. Host-adapted ecotype 
divergence was less than that recovered in pairwise popula- 
tion comparisons, especially when only significant compari- 
sons were considered. However, all significant population 
comparisons occurred between populations using alternate 
hosts suggesting some important host-population inter- 
action not obvious in the lingual mandible. Thus, in the 
occlusial mandible more reflective of structural features, re- 
sults suggest some influence of host plant adaptation, but 
only once the influence of other factors have been 
minimised (sex and rearing condition). These results suggest 
that host adaptation interacts extensively with population 
level differences and is specific to certain features accentu- 
ated in the different mandible orientations. 

Conclusions 

Collectively, presented data show that mandible shape in 
T. cristinae is under both genetic and environmental con- 
trol, mostly based on sexual differences, and to a lesser de- 
gree on differences between hosts and among populations. 
Mandible shape divergence in T, cristinae occurs predom- 
inantly along a different intraspecific factor than most other 
morphological traits and may, as a consequence, limit the 
progress toward ecological speciation between host-adapted 
ecotypes. Future work testing the causes of sexual dimorph- 
ism in trophic morphology in this species is required. 

Availability of supporting data 

The landmark data for each individual is stored through 
Labarchives.com and is avaflable at the foUowing link: 
http://dx.doi.org/10.6070/H4GQ6VPM. 
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Additional file 1: Supporting tests and comparisons. Supporting 
online information for this study outlines both common and global 
morphospace analyses, population level MU comparisons and 
relationships assessed between shape variation and the potential for 
gene flow among populations. 
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